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Abstract. Context: There is increasing need for good algorithms for modeling the aggregation and fragmentation of 
solid particles (dust grains, dust aggregates, boulders) in various astrophysical settings, including protoplanetary 
disks, planetary- and sub-stellar atmospheres and dense molecular cloud cores. Here we describe a new algorithm 
that combines advantages of various standard methods into one. 

Aims: The aim is to develop a method that 1) can solve for aggregation and fragmentation, 2) can easily include 
the effect and evolution of grain properties such as compactness, composition, etc., and 3) can be built as a 
coagulation/fragmentation module into a hydrodynamics simulation where it 3a) allows for non-'thermalized' 
non-local motions of particles (e.g. movement of particles in turbulent flows with stopping time larger than eddy 
turn-over time) and 3b) focuses computational effort there where most of the mass is. 

Methods: We develop a Monte-Carlo method in which we follow the "life" of a limited number of representative 
particles. Each of these particles is associated with a certain fraction of the total dust mass and thereby represents 
a large number of true particles which all are assumed to have the same properties as their representative particle. 
Under the assumption that the total number of true particles vastly exceeds the number of representative particles, 
the chance of a representative particle colliding with another representative particle is negligibly small, and we 
therefore ignore this possibility. This now makes it possible to employ a statistical approach to the evolution of 
the representative particles, which is the core of our Monte Carlo method. 

Results: The method reproduces the known analytic solutions of simplified coagulation kernels, and compares well 
to numerical results for Brownian motion using other methods. For reasonably well-behaved kernels it produces 
good results even for moderate number of swarms. 
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1. Introduction 

Dust particle aggregation is a very common process in 
various astrophysical settings. In protoplanetary disks the 
aggregation of dust particles forms the ve ry ini tial step 
of planet formation (see e.g. Dominik et al. l2007t h It also 
modifies the optical properties of the disk, and it has in- 
fluence on the che mistry and free electro n abu ndance in 
a disk (Sano et al. bOOOt Semenov et al. 12004 Ilgner & 
Nelson |200(J. The appearance and evolution of a proto- 
planetary disk is therefore critically affected by the dust 
aggregation process. In sub-stellar and planetary atmo- 
spheres the aggregation of dust particles and the coagula- 
tion of fluid droplets can affect the structure of cloud lay- 
ers. It can therefore strongly affect the spectrum of these 
objects and influence the local conditions within these at- 
mospheres. The process of aggregation/coagulation and 
the reverse process of fragmentation or cratering are there- 



fore important processes to understand, but at the same 
time they are extremely complex. 

Traditional methods solve the Smoluchowski equation 
for the particle mass distribution function f(m), where 
/(m) is defined such that f(m)dm denotes the number 
of particles per cubic centimeter with masses in the in- 
terval [m, m + dm]. This kind of method has been used in 
m any pa pers on dust coa gulation befo re (e.g. Nakagawaet 



al. 119811: Weiden schilli ng U98J; U993 Schmitt et al. 11997 : 
Suttner k, Yor ke Il999t Tanaka et al. 120051: Dullemond & 
Dominik [200ol Nomura & Nakagawa |200fj ). Methods of 
this kind are efficient, but have many known problems. 
First of all a coarse sampling of the particle mass leads 
to systematic error s such as the acceleration of growth 
(Ohtsuki et al. ll990h . High resolution is therefore required, 
which may make certain problems computationally expen- 
sive. Moreover, if one wishes to include additional prop- 
erties of a particle, such as porosity, charge, composition 
etc, then each of these properties adds another dimension 
to the problem. If each of these dimensions is sampled 
properly, this can quickly make the problem prohibitively 
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computationally expensive. Finally, the traditional meth- 
ods are less well suited for modeling stochastic behavior of 
particles unless this stochastic behavior can be treated in 
an averaged way. For instance, in protoplanetary disks if 
the stopping time of a particle is roughly equal to the tur- 
bulent eddy turn-over time, then the velocity of a particle 
with respect to the gas is stochastic: at the same loca- 
tion there can exist two particles with identical proper- 
ties but which happen to have different velocities because 
they entered the eddy from differ ent d irections (see e.g. 
the simulations by Johansen et al. 120061 ). 

To circumvent problems of this kind Ormel et 



al. (|2007[ ) have presented a Monte Carlo approach to coag- 
ulation. In this approach the particles are treated as com- 
putational particles in a volume which is representative 
of a much larger volume. The simulation follows the life 
of N particles as they collide and stick or fragment. The 
collision rates among these particles arc computed, and 
by use of random numbers it is then determined which 
particle collides with which. The outcome of the collision 
is then determined depending on the properties of the two 
colliding particles and their relative impact velocity. This 
method, under ideal conditions, provides the true simula- 
tion of the process, except that random numbers are used 
in combination with collision rates to determine the next 
collision event. This method has many advantages over 
the tradiational methods. It is nearly trivial to add any 
number of particle properties to each particle. There is less 
worry of systematic errors because it is so close to a true 
simulation of the system, and it is easy to implement. A 
disadvantage is that upon coagulation the number of com- 
putational particles goes down as the particles coagulate. 
Ormel et al. solve this problem by enlarging the volume 
of the simulation and hence add new particles, but this 
means that the method is not very well suited for model- 
ing coagulation within a spatially resolved setting such as 
a hydrodynamic simulation or a model of a protoplanetary 
disk. 

It is the purpose of the present paper to present an 
alternative Monte Carlo method which can quite nat- 
urally deal with extremely large numbers of particles, 
which keeps the number of computational particles con- 
stant throughout the simulation and which can be used in 
spatially resolved models. 

2. The method 

2.1. Fundamentals of the method 

The fundamental principle underlying the method we 
present here is to follow the behavior of a limited number 
of representative particles whose behavior is assumed to be 
a good representation of all particles. In this approach the 
number of physical particles N can be arbitrarily large. In 
fact it should be very much larger than the number of rep- 
resentative particles n, so that the chance that one repre- 
sentative particle collides with another representative par- 
ticle is negligible compared to the collisions between a rep- 



resentative particle and a non-representative particle. In 
other words, ii N ^> n, we only need to consider collisions 
between a representative particle and a non-representative 
particle. The number of collisions among representative 
particles is too small to be significant, and the collisions 
among non-representative particles are not considered be- 
cause we focus only on the behavior of the representative 
particles. 

Suppose we have a cloud of dust with N = 10 20 physi- 
cal particles, with a specific size distr ibutio n, for instance, 
MRN (Mathis, Rumpl & Nordsieck Il977t ). Let the total 
mass of all these particles together be M to t and the vol- 
ume be V. We randomly pick n particles out of this pool, 
where n is a number that can be handled by a computer, 
for instance, n — 1000. Each representative particle i has 
its own mass m, and possibly other properties such as 
porosity pi or charge Ci assigned to it. We now follow the 
life of each of these n = 1000 particles. To know if repre- 
sentative particle i = 20 collides with some other object, 
we need to know the distribution function of all physical 
particles with which it can collide. However, in the com- 
puter we only have information about the n representative 
particles. We therefore have to make the assumption that 
the distribution function set up by the n representative 
particles is representative of that of the N physical parti- 
cles. We therefore assume that there exist 

M tot 

n k = 77 (1) 

nm k V 

physical particles per cubic cm with mass m k , porosity p k , 
charge Ck etc., and the same for each value of k, including 
k = i. In this way, by assumption, we know the distribu- 
tion of the N physical particles from our limited set of n 
representative particles. One could say that each represen- 
tative particle represents a swarm of M io t/nm,i physical 
particles with identical properties as the representative 
one. One could also say that the true distribution of N 
particles is, by assumption, that of the n representative 
ones. The rate of collisions that representative particle i 
has with a physical particle with mass m k etc. is then: 

M tot 



n(k) = n k a ik Av ik 



ik 



(2) 



nrrikV 

where Ci k is the cross-section for the collision between 
particles with properties i and fc, and Av ik is the average 
relative velocity between these particles. The total rate 
of collisions that representative particle i has with any 
particle is then: 



k 



r-i(fc) 



(3) 



and the total rate of collisions of any representative par- 
ticle is 

(4) 



The time-evolution of the system is now done as fol- 
lows. Let to be the current time. We now randomly choose 
a time step St according to: 

5t = — log(ran(seed)) (5) 
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where ran(seed) is a random number uniformly distributed 
between and 1. This means that a collision event hap- 
pens with one of the representative particles at time 
t = to + St. The chance P(i) that the event happens to 
representative particle i is: 



P(i) 



n 

r 



(6) 



So we can choose, using again a random number, which 
representative particle i has undergone the collision event. 
We now need to determine with what kind of physical par- 
ticle it has collided. Since the distribution of physical par- 
ticles mirrors that of the representative ones, we can write 
that the chance this particle has collided with a physical 
particle with properties k is: 



P(k\i) 



n{k) 



(7) 



With another random number we can thus determine 
which k is involved in the collision. Note that k can be 
i as well, i.e. the representative particle can collide with 
a physical particle with the same properties, or in other 
words: a representative particle can collide with a particle 
of its own swarm of physical particles. 

Now that we know what kind of collision has hap- 
pened, we need to determine the outcome of the collision. 
The most fundamental part of our algorithm is the fact 
that only representative particle i will change its proper- 
ties in this collision. Physical particle k would in principle 
also do so (or in fact becomes part of the new represen- 
tative particle), but since we do not follow the evolution 
of the physical particles, the collision will only modify the 
properties of representative particle i. By assumption this 
will then automatically also change the properties of all 
physical particles associated with representative particle i. 
Statistically, the fact that the particles k are not modified 
is "corrected for" by the fact that at some point later the 
representative particle k will have a collision with physical 
particle i, in which case the properties of the k particles 
will be modified and not those of i. This then (at least in 
a statistical sense) restores the "symmetry" of the inter- 
actions between i and k. If the collision leads to sticking, 
then the resulting particle will have mass m = rrii + . 
This means that representative particle i will from now 
on have mass <— m, + m^. Representative particle 
k is left unaffected as it is not involved in the collision. 
The interesting thing is now that, because by assump- 
tion the representative particle distribution mirrors the 
real particle distribution, the swarm of physical particles 
belonging to the modified representative particle i now 
contains fewer physical particles, because the total dust 
mass M = Mt t/n of the swarm remains constant. 

If a collision results in particle fragmentation, then the 
outcome of the collision is a distribution function of debris 
particles. This distribution function can be written as a 
function fd(m) of debris particle mass, such that 



and the function fd,(m) has to be determined by labora- 
tory experiments or detailed computer simula tions of in- 
dividual particle collisions (see Dominik et al. 120071 for a 
review). The new value of rrii for the representative parti- 
cle is now randomly chosen according to this distribution 
function by solving the equation 



mfd{m)dm — ran(seed)(mj + rrif.) 



(9) 



for fh and assigning m, <— m. In other words: we randomly 
choose a particle mass from the debris mass distribution 
function, i.e. the choice is weighed by fragment mass, not 
by fragment particle number. This can be understood by 
assuming that the true representative particle before the 
collision is in fact just a monomer inside a larger aggre- 
gate. When this aggregate breaks apart into for instance 
one big and one small fragment it is more likely that this 
representative monomer resides in the bigger chunk than 
in the smaller one. 

After a fragmenting collision the rrii will generally be 
smaller than before the collision. This means that the 
number of physical particles belonging to representative 
particle i increases accordingly. Note that although the 
collision has perhaps produced millions of debris particles 
out of two colliding objects, our method only picks one 
of these debris particles as the new representative particle 
and forgets all the rest. Clearly if only one such destructive 
collision happens, the representative particle is not a good 
representation of this entire cloud of debris products. But 
if hundreds such collisions happen, and are treated in the 
way described here, then the statistical nature of Eq. ((9]) 
ensures that the debris products are well represented by 
the representative particles. 

The relative velocity Av can be taken to be the aver- 
age relative velocity in case of random motions, or a sys- 
tematic relative velocity in case of systematic drift. For 
instance, for Brownian motion there will be an average 
relative velocity depending on the masses of both par- 
ticles, but differential sedimentation in a protoplanetary 
disk or planetary atmosphere generates a systematic rela- 
tive velocity. Also, for the Brownian motion or turbulent 
relative velocity one can, instead of using an average rel- 
ative velocity, choose randomly from the full distribution 
of possible relative velocities if this is known. This would 
allow a consistent treatment of fluctuations of the relative 
velocities which could under some circumstances become 



important (see e.g. Kostinski & Shaw l2005T ). 



mfd{m)dm — rrii + TTik 



(8) 



2.2. Computer implementation of the method 

We implemented this method in the following way. For 
each of the representative particles we store the mass m, 
and all other properties such as porosity, charge, compo- 
sition etc. Before the start of the Monte Carlo procedure 
we compute the full collision rate matrix and we 

compute the r*j as well as r. For these collision pairs (i, k) 
we now have to determine the cross section of particles 
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as well as their systematic relative velocity, such as differ- 
ent drift speeds, and the random relative velocity, such as 
Brownian or turbulent motion. The random motions can 
be determined with a random number from the relative 
velocity probability distribution function if that is known. 
If that is not known in sufficient detail, one can also take 
it to be the average relative velocity, for which more often 
analytic formulae exist in the literature. 

We determine beforehand at which times t sav .„ we 
want to write the resulting rrii and other parameters to 
a file. The simulation is now done in a subroutine with a 
do-while loop. We then determine St using a random num- 
ber (see Eq. [5]), and check if t + 5t < t sav ^ n , where t saVtTl 
is the next time when the results will have to be stored. 
If t + St < i S av,nj then a collision event occured before 
t S a.v,n- We will handle this event according to a procedure 
described below, we set t <^ t + St and then return to 
the point where a new St is randomly determined. If, on 
the other hand, t + St > i saVjn then we stop the proce- 
dure, return to the main program and set t <— t S av,n- The 
main program can then write data to file and re-call the 
subroutine to a time i S av,n+i or stop the simulation alto- 
gether. Note that when the subroutine is called again for 
a next time interval, it does not need to know the time of 
the previously randomly determined event which exceeded 
isav,n- Of course, one could memorize this time and take 
that time as the time of the next event in the next time 
interval, but since the events follow a Poisson distribution, 
we do not need to know what happened before t sav .„ to 
randomly determine the new time t + St of the next event. 

Now let us turn to what happens if a collision event oc- 
curs, i.e. occurs between time t and t+St. We then first de- 
termine which representative particle i is hit, which is done 
by generating a random number and choosing from the 
probability distribution of collision rates, as described in 
Section r2.ll Similarly we determine the non-representative 
particle with which it collides, or in other words: we de- 
termine the index k of the "swarm" in which this non- 
representative particle resides. Finally, we must determine 
the impact parameter of the collision, or assume some av- 
erage impact parameter. 

Now we employ a model for the outcome of the colli- 
sion. This is the collision subroutine of our Monte Carlo 
method. It is here where the results of laboratory exper- 
iments come in, and the translation of such experiments 
into a coagulation kernel is a major challenge which we do 
not cover here. The collision model must be a quick for- 
mula or subroutine that roughly represents the outcomes 
of the detailed laboratory collision experiments or detailed 
numerical collision models. It will give a probability func- 
tion fd for the outcoming particle masses and properties. 
From this distribution function we pick one particle, and 
from this point on our representative particle i will attain 
this mass and these properties. The collision partner k will 
not change, because it is a non-represctativc particle from 
that swarm that was involved in the collision, and we do 
not follow the life of the non-representative particles. We 
therefore ignore any changes to that particle. 



We now must update Vj(t) for all I with fixed j = i 
and for all j with fixed I = k: we update a row and a 
column in the rj(l) matrix. Having done this, we must also 
update Tj for all j. This would be an n 2 process, which is 
slow. But in updating r.j (I) we know the difference between 
the previous and the new value, and we can simply add 
this difference to Tj for each j. Only for j = i we must 
recompute the full rj again, because there all elements 
of that row have been modified. Using this procedure we 
assure that we limit the computational effort to only the 
required updates. 

2.3. Acceleration of the algorithm for wide size 
distributions 

One of the main drawbacks of the basic algorithm de- 
scribed above is that it can be very slow for wide size dis- 
tributions. Consider a swarm of micron sized dust grains 
that are motionless and hence do not coagulate among 
each other. Then a swarm of meter sized boulders moves 
through the dust swarm at a given speed, sweeping up the 
dust. Let us assume that also the boulders are not collid- 
ing among each other. The only mode of growth is the 
meter-sized boulders sweeping up the micron sized dust. 
For the boulder to grow a factor of 2 in mass it will have 
to sweep up 10 18 micron sized dust particles. Each impact 
is important for the growth of the boulder, but one needs 
10 18 such hits to grow the boulder a factor of 2 in mass. 
The problem with the basic algorithm described above is 
that it is forced to explicitly model each one of these 10 18 
impacts. This is obviously prohibitively expensive. 

The solution to this problem lies in grouping collisions 
into one. Each impact of a dust grain on a boulder only in- 
creases the boulder mass by a minuscule fraction. For the 
growth of the boulder it would also be fine to lower the 
chance of an impact by 10 16 , but (fit happens, then 10 16 
particles impact onto the boulder at once. Statistically this 
should give the same growth curve, and it accelerates the 
method by a huge factor. However, it introduces a fine- 
tuning parameter. We must specify the minimum increase 
of mass for coagulation (dm max ). If we set dm max to, for 
instance, 10%, then we may expect that the outcome also 
has errors of the order of 10%. This error arises because 
by increasing the mass of the bigger body in steps of 10%, 
we ignore the fact that the mass at some time should in 
fact be somewhere in between, which cannot be resolved 
with this method. This is, however, not a cumulative er- 
ror. While the mass of the bigger body may sometimes 
be too low compared to the real one, it equally probably 
can be too large. On average, by the Poisson nature of the 
collision events, this averages out. But it is clear that the 
smaller we take this number, the more accurate it becomes 
- but also the slower the method becomes. It is therefore 
always a delicate matter to choose this parameter, but for 
problems with a large width of the size distribution this 
acceleration is of vital importance for the usability of the 
method. 
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2.4. Including additional particle properties 

We mentioned briefly the possibility of adding more parti- 
cle properties to each representative particle. This is very 
easy to do, and it is one of the main advantages of a Monte 
Carlo method over methods that directly solve the integral 
equations for coagulation. One of the main properties of 
interest to planet formation is porosity or fractal structure 
of the aggregate. Two aggregates with the same mass can 
have vastly different behavior upon a collision if they have 
different compactness. A fluffy aggregate may break apart 
already at low impact velocities while a compact aggre- 
gate may simply bounce. Upon collisio ns the se properties 
may in fact also change. Ormel et al. (|20071 ) studied the 
effect of porosity and how it changes over time, and they 
also used a Monte Carlo approach for it. 

If one wishes to include particle properties in a tra- 
ditional method which solves the integral equations of 
coagulation (the Smoluchowski equation), then one in- 
creases the dimensionality of the problem by 1 for each 
property one adds. With only particle mass one has a 
distribution function f(m,t) while adding two particle 
properties p\ and pi means we get a distribution func- 
tion f(m,pi,p2,i), making it a 4-dimensional problem. 
Methods of this kind must treat the complete phase 
space spanned by (m,pi,p2,t). This is of course possi- 
ble, but c omput ationally it is a very challenging task (see 
Ossenkopf ll993f ). In contrast, a Monte Carlo method only 
sparsely samples phase space, and it samples it only there 
where a significant portion of the total dust mass is. A 
Monte Carlo method focuses its computational effort au- 
tomatically there where the action is. The drawback is 
that if one is interested in knowing the distribution func- 
tion there where only very little mass resides, then the 
method is inaccurate. For instance, in a protoplanetary 
disk it could very well be that most of the dust mass is 
locked up in big bodies (larger than 1 meter) which are 
not observable, and only a promille of the dust is in small 
grains, but these small grains determine the infrared ap- 
pearance of the disk because they have most of the solid 
surface area and hence most of the opacity. In such a case 
a Monte Carlo method, by focusing on where most of the 
mass is, will have a very bad statistics for those dust grains 
that determine the appearance of the disk. For such goals 
it is better to use the traditional methods. But if we are in- 
terested in following the evolution of the dominant portion 
of the dust, then Monte Carlo methods naturally focus on 
the interesting parts of phase space. 

3. Discussion of the method 

3.1. Conservation of particle number 

There are a few peculiarities of the method described here 
that may, at first sight, appear inconsistent, but are sta- 
tistically correct. For instance if we return to our exam- 
ple of a swarm of tiny particles and a swarm of boulders, 
i.e. n = 2 with representative particle 1 being a micron 
sized particle and representative particle 2 being a meter 



sized particle, then we encounter an apparent paradox. We 
again assume that collisions only take place between 1 and 
2, but not between 1 and 1 or 2 and 2. The chance that 
representative particle 1 hits a meter size particle is much 
smaller than the chance that representative particle 2 hits 
a micron size particle. What will happen is that represen- 
tative particle 2 will have very many collision events with 
small micron size grains, and thereby slowly and gradu- 
ally grows bigger, while representative particle 1 will only 
have a collision with big particle after a quite long time 
and immediately jumps to that big size. While representa- 
tive particle 2 grows in mass, the number of big physical 
particles decreases in order to conserve mass. This may 
seem wrong, because in reality the number of big boul- 
ders stays constant, and these boulders simply grow by 
sweeping up the small dust. The solution to this paradox 
is that the average time before representative particle 1 
hits a big (k = 2) particle is of the same order as the time 
it takes for representative particle 2 to grow to twice its 
mass by collecting small particles. So, very roughly, by the 
time the big particle has doubled its mass, and therefore 
the number of physical particles belonging to k = 2 has 
reduced by 50%, the representative particle 1 has turned 
into a big particle, corresponding, statistically, to the other 
50% of big particles that was missing. If we are a bit more 
precise, the statistics do not add up precisely in this way if 
we have only 1 swarm of small and 1 swarm of big bodies. 
If, however, one has N swarms of small and N swarms of 
big bodies, and again assume that only the big bodies can 
sweep up the smaller ones, then if N 3> 1 the statistics 
adds up perfectly: one finds that after all the growth has 
taken place, the average mass of the bodies is twice that 
of the original big bodies. In Figure Q] we do precisely this 
experiment, and the left panel shows that for N — 500 the 
mass distribution of the big bodies averages to the right 
value, albeit with a spread of 10% FWHM while in reality 
this spread should be 0. The right panel shows how the 
average final mass depends on N. For small N the statis- 
tics clearly do not add up, but for large N they do and 
produce the right value (final mass is twice initial mass of 
the big bodies). So statistically the number of big particles 
is restored to the correct value, but there is then unfor- 
tunately still a large statistical noise on it. The particle 
number is therefore not exactly conserved in our method, 
but statistically it is. 



3.2. The number of representative particles 

It is obvious that for high number of representative parti- 
cles n we will get better results than for low n. But there 
are two issues here. First of all, the higher n, the better the 
representative particles represent the true physical distri- 
bution of particles. For problems that result in wide size 
distributions this is all the more crucial. An inaccurate 
representation of the true size distribution could lead to 
systematic errors. But another reason for taking a high n 
is simply because we want our end-result to have as little 
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Fig. 1. Results of the test problem with N swarms of small particles and N swarms of big bodies, as discussed in 
Section I37T1 Left: histogram of the final masses of the bodies relative to the initial mass of the big bodies, for N = 500. 
Right: average mass relative to the initial mass of the big bodies as a function of N. 



as possible noise. If the result is too noisy, then it is useless. 
Taking n too big, however, makes the code slow because 
more representative particles have to be followed, and for 
each of these particles we must check for a larger number 
of possible collision partners k. The problem scales there- 
fore as n 2 . If the expected size distributions are not too 
wide, one can use an intermediately large number of repre- 
sentative particles, say n, for the simulation, but redo the 
simulation m times such that n = mri, and average the 
results of all m si mulat ions. This approach was also used 
by Ormel et al. (|2007l ). This gives the same amount of 
noise on the end-result, but scales as n 2 m — n 2 /m, which 
is m times faster than the n 2 scaling. This works, how- 
ever, only if the coagulation/fragmentation kernel is not 
too sensitive to the exact distribution of collision partners. 

Interestingly, if the kernel is very insensitive to the ex- 
act distribution of collision partners, then, in principle, 
one could run the model with only a single representative 
particle n = 1, because the collision partner of represen- 
tative particle i could be equal to k = i. Of course, a 
single representative particle means that we assume that 
all physical particles have the same size, or in other words: 
that we have an infinitely narrow size distribution. 

To decide about the sufficient number of representa- 
tive particles, one has to compare the results of the MC 
code with the analytical solutions of the three test ker- 
nels (see Section 01 . In a given time of the simulation the 
mean mass and the shape of the distribution function for 
all three test kernels must be followed accurately. It is es- 
pecially important to reproduce the linear and product 
kernels accurately as the realistic kernels of dust particles 
are similar to these. 

Of course the progression from 'not sufficient' and 'suf- 
ficient' number of representative particles is smooth and 
in general the more representative particles we use, the 
more accurate the produced result will be. The sufficient 
number of representative particles (n) as given in Section 



U) are only suggestions, the error of the distribution func- 
tions were not quantified. 



3.3. Limitations of the method 

One of the fundamental limitations of the method de- 
scribed here is that we assume N ^> n. We can model 
the growth of particles by coagulation in a protoplanetary 
disk or in a cloud in a planetary atmosphere, but we can 
not follow the growth to the point where individual large 
bodies start to dominate their surroundings. For instance, 
if we wish to follow the growth of dust in a protoplanetary 
disk all the way to small planets, then the method breaks 
down, because N is then no longer much bigger than n, 
and interactions among representative particles become 
likely. Also, for the same reason, run-away grow th prob- 
lems such as electrostatic gelation (Mokler 12007 ) cannot 
be modeled with this method. 



Another limitation is encountered when modeling 
problems with strong growth and fragmentation happen- 
ing at the same time. This leads to very wide size distri- 
butions, and the typical interval between events is then 
dominated by the smallest particles, whereas we may be 
interested primarily in the growth of the biggest particles. 
In such a situation a semi-steady-state can be reached 
in which particles coagulate and fragment thousands of 
times over the life time of a disk. The Monte Carlo method 
has to follow each of these thousands of cycles of growth 
and destruction, which makes the problem very "stiff" . 
Methods using the integral form of the equations, i.e. 
the Smoluchowski equation, can be programmed using 
implicit integration in time so that time steps can be 
taken which are much larger than the typical time scale of 
one growth- fr agmen tation cycle without loss of accuracy 
(Brauer et al. l2008h . This is not possible with a Monte 
Carlo method. 
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Fig. 2. Test against the constant kernel [Kij = 
1). The particles were binned and the distribution 
function was produced at dimensionless times t — 
0, 10°, 10\ 10 2 , 10 3 , 10 4 , 10 5 . The dashed lines show the 
analytical solution. This run was produced by simulat- 
ing 200 representative particles five times and producing 
the average of these. In this case dm max is 0.1. 



Fig. 3. Test against the linear kernel (Kij = m,i + rrij). 
The particles were binned and the distribution function 
was produced at dimensionless times t = 0, 4, 8, 12, 16, 20. 
The dashed lines show the analytical solution. This run 
was also produced by simulating 200 representative parti- 
cles five times and producing the average of these. In this 
case dm max is 0.1. 



4. Standard tests and results 

In this section we test our coagulation model with kernels 
that have analytical solutions. Furthermore we show the 
first results of applying this model to protoplanetary disks 
introducing Brownian motion and turbulence induced rel- 
ative velocities as well as a new property of dust particles 
n amely the porosity (or enlargement factor, see Ormel et 
al 120071 ). and a simple fragmentation model. 

To follow dust coagulation and fragmentation, one has 
to follow the time evolution of the particle distribution 
function at a given location in the disk (/(y,t)), where 
y contains the modeled properties of the dust grains, in 
our case these will be the mass (to) and the enlargement 
factor (*), f(m, \P, t). 

In most of the coagulation models so far the only used 
dust-property was the particle mas s. Then one can use th e 
so called Smoluchowski equation ( Smoluchowskil ( 19161) ) 
to describe the time-evolution of /(to): 

= -/(m) J dm'K(m,m')f(m') + 

i J dm! Kim', to - m')f(m')f(m - to'). (10) 

The first term on the right hand side represents the loss 
of dust in the mass bin m by coagulation of a particle 
of mass to with a particle of mass to'. The second term 
represents the gain of dust matter in the mass bin to by 
coagulation of two grains of mass to' and to — to'. K is 
the coagulation kernel, it can be written as 

K(mi,m 2 ) = er c (mi,TO 2 ) x Av(m 1 ,m 2 ), (11) 

the product of the the cross-section of two particles and 
their relative velocity. We consider all the three kernels for 



which there exist analytical solutions: The constant kernel 



(Kij — 1), the linear kernel (Kij 



j) and the 



product kernel (Kij 
are de scribed e.g. in 
(jl990h . 



= mi x to,). The anal ytical solutions 
Ohtsuki et all (|l990h and Wetherill 



We test our method against these three kernels, leaving 
the enlargement factor unchanged, always unity. Further 
important properties of the dust particles, such as ma- 
terial density and volume density, are also always unity. 
The (dimensionless) time evolution of the swarms is fol- 
lowed and at given times the particles are binned by mass 
so that we can produce /(m). On Figures [2] and [3] the 
y axis shows /(to) x to 2 , the mass density per bin. The 
analytical soluti ons, taken from I Ohtsuki et all (|l990l) and 
Wetherill (|l990f) . are overplotted with dashed line. The 
number of particles were chosen to be n = 200, to = 5, so 
altogether 1000 representative particles were used in the 
model except for the product kernel where more represen- 
tative particles were used to achieve better results. 

In the case of the constant kernel (Figure[2]), we started 
our simulation with MRN size distribution (n(a) oc a -3 ' 5 ), 
the results were saved at t = 0, 10°, 10 1 , 10 2 , 10 3 , 10 4 , 10 5 . 
It is interesting to note that this kernel is not sensitive 
to the initial size distribution. As the system evolves, it 
forgets the initial conditions. Another interesting property 
of this kernel that our model can reproduce the analytical 
solution even with very limited number of representative 
particles (even for n = 5!) but of course with higher noise. 
It is possible to use only one representative particle, which 
means that the representative particle collides with par- 
ticles from its own swarm which basically results in pure 
CCA growth (Cluster-Cluster Aggregation). Interestingly, 
the mean mass of the distribution function is followed cor- 
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10000 
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proximately n « 500 representative particles to correctly 
reproduce it. 

The required CPU time for these test cases is very low, 
some seconds only. 

We conclude that our Monte Carlo method reproduces 
the constant and linear test kernels without any problem 
even with low number of representative particles. On the 
other hand the method has difficulties with the product 
kernel, but before the formation of the 'run-away' particle, 
we can reproduce the kernel. The relatively low number 
of representative particles needed to sufficiently reproduce 
the test kernels is very important for future applications 
where whole disk simulations will be done and there will 
likely be regions containing low numbers of particles. 



Fig. 4. Test against the product kernel (Kij = rrii x irij). 
The particles were binned and the distribution function 
was produced at dimensionless times t = 0.4, 0.7, 0.95. The 
dashed lines show the analytical solution. This run was 
produced by simulating 1000 representative particles ten 
times and producing the average of these. In this case 
dm max is set to be 0.05. 



rectly but the shape of the function changes, additional 
spikes appear on it. 

The linear kernel is known to be more problematic be- 
cause the mean mass of the particles grows exponentially 
with time. Our model, however reproduces this kernel very 
well, too, as it can be seen in Figure [3] The results were 
saved at t = 0, 4, 8, 12, 16, 20. We note that using low num- 
ber of representative particles with this kernel also works 
relatively well, the minimum number of swarms needed 
to reproduce the exponential time evolution of the mean 
mass is h w 100. This is larger than for the constant ker- 
nel. It shows that for the linear kernel collisions between 
particles of unequal mass are contributing significantly to 
the growth, whereas for the constant kernel the growth is 
dominated by collisions between roughly equal size parti- 
cles. Using n <C 100 results in distorted distribution func- 
tion: neither the mean mass nor the actual shape of the 
distribution function is correct. 

The product kernel is the hardest to reproduce. The 
peculiarity of this kernel is the following: Using dimen- 
sionless units, a 'run-away' particle is produced around 
t = 1 , which collects all the o ther particles present in 
the simulation (Wetherill Il990h . The difficulty arises m 
our Monte Carlo code when the mass of the representa- 
tive 'run-away' particle reaches the mass of its swarm. In 
other words, the number of physical particles belonging to 
the representative 'run-away' particle is close to unity. In 
this case the original assumption of our method (we only 
need to consider collisions between a representative parti- 
cle and a physical particle) is not valid anymore. However, 
as Figure |4] shows, we can relatively well reproduce this 
kernel before t = 1. In the case of this kernel, we need ap- 



5. Applications to protoplanetary disks 

We use the Monte Carlo code to follow the coagulation 
and fragmentation of dust particles in the midplanc of a 
protoplanetary disk at 1 AU from the cent ral star. Our 
disk m odel is identical with the one used bv iBrauer et al 



( 20071 ). We proceed step by step. First relative velocities 
induced by Brownian motion and turbulence without the 
effects of porosity are included fSec. 15. 1| ). 

The next step is to include a fragmentation model 
(Sec. EH). 

In the final step porosity is included (Sec. I5.3[) . W e 
use the porosity model described in Ormel et al. ~ i|2007l) . 
At this point we co mpare and check again our code with 
Ormel et al. (2007]) using their input parameters but not 
including the rain out of particles. 

5.1. Relative velocities 

We include two processes in calculating the relative veloc- 
ities: Brownian motion and turbulence. 

Brownian motion strongly depends on the mass of the 
two colliding particles. The smaller their masses are, the 
more they can be influenced by the random collisions with 
the gas molecules/atoms. One can calculate an average 
velocity given by 



Av B {mi,m 2 ) 



t 8kT(m 1 + m 2 ) 
7rmi m2 



(12) 



For micron sized particles, relative velocity can be in the 
order of magnitude of 1 cm/s, but for cm sized particles 
this value drops to 10 -7 cm/s. If growth is only governed 
by Brownian motion, it leads to very slow coagulation, a 
narrow size distribution and fluffy dust particles, so called 
cluster-cluster aggregates (CCA). 

The gas in the circumstellar disk is turbulent, thus the 
dust particles experience acceleration from eddies with 
different sizes and turnover tim es. Th is process is very 
complex, but Ormel and Cuzzi (|2007Q provided limiting 
closed-form expressions for average relative turbulent ve- 
locities between two dust particles. Their results are also 
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The fragmentation energy is then defined as follows: 



-o 4 



Particle radius cm 



Fig. 5. The relative velocity caused by Brownian motion 
and turbulence for different sized particles. The black line 
shows the fragmentation barrier. Collision events situated 
between these two lines result in fragmentation if porosity 
is not included. Physical parameters of the disk: the dis- 
tance from the central star is 1 AU, temperature is 200 K, 
the density of the gas is 8.73 x 10 12 particle/cm 3 , and the 
turbulent parameter, a — 10~ 3 . Parameters of the dust: 
monomer radius is ao = 0.4/im, material density is p — 1.6 
g/cm 3 . 

valid for particles with high Stokes numbers. They dis- 
tinguished three regimes: a.) the stopping times of both 
dust particles are smaller than the smallest eddy-turnover 
time {t\,t2 < t v , tightly coupled particles); b.) the stop- 
ping time is between the smallest and largest turnover 
time (t v < t\ < ti,, intermediate regime); c.) the stopping 
time is bigger than the largest turnover time (t\ > tj, , 
heavy particles). For details see Ormel and Cuzzi (|2007l ). 
We used a — 10~ 3 for the turbulence parameter. 

To illustrate the relative velocity of dust particles with- 
out the effects of porosity, we provide Figure [5] This con- 
tour plot includes Brownian motion and turbulent relative 
velocities. The Brownian motion is negligible for particles 
bigger than 10~ 2 cm. 

5.2. Fragmentation model 

The collision energy of the particles is 



E=~ '—^Av 2 = -pAv 2 , 

2 mi + mi I 



(13) 



where p is the reduced mass. We need to define some quan- 
tities of the dust particles. E ro u is the rolling energy of two 
monomers. For mon omers of the same siz e it is given by 
(Dominik & Tielens |l997| Blum & Wurm|2000|) 



E 



roll 



-7raoF ro u, 



(14) 



where oo is the monomer radiu s, F ro u is the rolling force 
measured by Heim et al. (1999). Its value is F roU = (8.5 ± 
1.6) x 10~ 5 dyn for Si0 2 spheres. 



Efrag — N c X Ebreak — 3iV X E ro l, 



(15) 



where N c is the total number of contact surfaces between 
monomers (for simplicity it is taken to be 3N, where N 
is the number of monomers in the particle), Eb rea k is the 
energy needed to break the bond between two monomers 
(its order of magnitude is similar to E ro u for these param- 
eters). 

If the collision energy of two particles is higher than 
the corresponding fragmentation energy, then the aggre- 
gate is destroyed and monomers are produced. Note that 
although assuming a complete destruction of the collided 
dust particles, we are interested in the critical energy 
where the first fragmentation event happens. This is the 
reason why the fragmentation energy is assumed to be 
lower than the energy needed for catastrophic fragmenta- 
tion. It is a simplification of the model to assume that the 
debris particles will be monomers. This is a very simpli- 
fied fragmen tation model used previously by Dullemond 
Dominik (2005). A more realistic model would be the 



one used bv lBrauer et al. I (|2007l) . 

We show the fragmentation barrier in Figure [5] with 
black lines. If collision happens in the regime between 
these two lines, that results in fragmentation. 

5.2.1. Results 

A simulation was made including these effects in a spec- 
ified location of the disk. We choose the location to be 
1 AU distance fr om the centra l solar type star. Using 
the disk model of iBrauer et al.l ( 2007 ). the temperature 
at this distance is approximately 200 K, the density of the 
gas is 8.73 x 10 12 cm~ 3 , the gas-to-dust ratio is 100 and 
we choose the turbulent parameter to be a = 10~ 3 , the 
Reyn olds number is Re = 10 8 (based on Ormel & Cuzzi 
20071 ). The dust monomers have the following properties: 
the monomer radius is ao = 0.4 /im, material density is 
p = 1.6 g/cm 3 . With the used parameters the fragmen- 
tation velocity is Avf rag ~ 8 m/s, though it is some- 
what larger for equal sized agglomerates. It is important 
to note that this value is very sensitive to the monomer ra- 
dius (ao) and material density (ft), because smaller /lighter 
monomers mean more contact surfaces (higher N for the 
same mass) and therefore higher fragmentation energy. 

Using these input parameters we simulated the evolu- 
tion of the dust particles for 3 x 10 3 years so that we reach 
an equilibrium between coagulation and fragmentation. 
Figure [5] shows the resulting normalized size distributions 
in times after t = 3 x 10°, 3 x 10 1 , 3 x 10 2 and 3 x 10 3 years. 
We used n = 100 particles averaging over m — 100 times 
(10 4 particles altogether). The required CPU time to per- 
form this simulation is 1.5 hours approximately. dm max is 
set to be 0.001 from now on in every simulation. We would 
like to note that giving dm max fSection l2.3[) a higher value 
would decrease the CPU time. 

One can see that coagulation happens due to Brownian 
motion in the beginning of the simulation (until 3 x 10 1 
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Fig. 6. The evolution of dust particles including the ef- 
fects of Brownian motion and turbulence. Porosity is not 
included in this model. The particle distribution is saved 
after t = 3 x 10° years - dash-dot line, 3 x 10 1 



years 



dashed line, 3 x 10 years - dotted line, and 3 x 10 years 
- continuous line. 



years) but after that turbulence takes over and the first 
fragmentation event happens after roughly 10 3 years. 
After this event the "recycled" monomers start to grow 
again, but as we see in Figure [5l particles can not reach 
bigger sizes than 0.07 cm. 

We would like to draw attention to the sudden de- 
crease of particles around 0.002 cm in Figure El This is 
the result of the turbulent relative velocity model used 
here (discussed in Sect. 15. 1|) . At this point the particles 
leave the 'tightly coupled particles' regime and enter the 
'intermediate' regime. But the transition in relative veloc- 
ity between these regimes is not smooth, there is a jump 
in relative velocity from ~20 cm/s to ^60 cm/s. As a 
result, particles coagulate suddenly faster and leave this 
part of the size distribution rapidly. Similar 'valleys' can 
be seen in the following figures with porosity, but the fea- 
ture is less distinct as the stopping times can be different 
for particles with same mass. 

5.3. Porosity 

To be able to quantitatively discuss the effect of porosity, 
we have to define the enlarg ement parameter following 
the discussion of Ormel et al. (|2007 ). If V is the extended 
volume of the grain and V* is the compact volume, than 
one can define the enlargement parameter ('J) as 



V 



(16) 



Compact volume is the volume occupied by the monomers 
not taking into account the free space between the 
monomer spheres. One can think of it as melting all the 
monomers into a single sphere, the volume of this sphere 
is the compact volume. We use compact radius later on, 



Fig. 7. The evolution of dust particles including the effects 
of Brownian motion and turbulence. Porosity is included 
in this model! The x axis shows the compact radius. The 
particle distribution is saved after t = 3 x 10° years - dash- 
dot line, 3 x 10 1 years - dashed line, 3 x 10 2 years - dotted 
line, and 3 x 10 3 years - continuous line. Note that the 
scaling of the x axis is different from Figure El 



which is the radius of this sphere. In the previous sec- 
tion the mass/volume ratio was constant for the particles. 
Therefore we could automatically calculate the mass of 
the particle if the radius was known or vice-versa. But 
from now on a particle with given mass m can have a 
wide range of effective radii depending on its enlargement 
parameter. 

It is essential to know how the enlargement parameter 
changes upon collisions. We have to refine our fragmenta- 
tion model and introduce two more regimes regard ing to 
collision energy. We use the model of Ormel et al. (|2007r ) 
and we only summarize their model here. 

The first regime is the low collision energy regime, 
where the collision energy is smaller than the restructuring 
energy (E < E restr , where E restr = 5E ro u), meaning that 
the particles stick where they meet, the internal structure 
of the grain does not change. 

The recipe for the resulting enlargement factor after 
the collision of two particles assuming that mi > ni2 then 
is 

V mi^xj 

where ($f) m is the mass averaged enlargement factor of 
the colliding particles: 



(*)r 



mi$i + TO2 v I / 2 

mi + m 2 



(18) 



Furthermore Scca is the CCA-characteristic exponent 
calculated b y det ailed numerical studies such as Paszun 
& Dominik (2006) {Scca — 0.95). ^ a dd is a necessary ad- 
ditional factor for the enlargement factor (for details see 
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Ormel et al.|2007): 



dd 



mi 
mi 



4"i exp 



10mo 



(19) 



where mo is the monomer mass. 

The second regime is the regime of compaction. The 
internal structure of the monomers inside the particle 
changes, this causes a decreasing porous volume. If the 
collision energy E res t r < E < Ef rag , we talk about com- 
paction. In this case the porosity after the collision be- 
comes 

* = (l-/c)((*) TO -l) + l, (20) 

where fc = E/(NE ro ii) — —AV/V is the relative com- 
paction. One can see that fc has to be smaller than unity 
otherwise $ in Eq. l20l becomes less than unity. But it can 
theoretically happen that E > NE ro u . In this long 
as the total collision energy remains below the fragmen- 
tation threshold, we assume that after compaction this 
excess energy goes back into the kinetic energy of the two 
colliding aggregates. The two aggregates therefore com- 
pactify and bounce, without exchanging mass or being 
destroyed. Bouncing is therefore included in this model, 
albeit in a crude way. 

The third regime is fragmentation as it was discussed 
in the previous section (E > Ef rag ). We use the same frag- 
mentation model as before so the result of a fragmenting 
collision arc monomers. 
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Fig. 8. This figure shows the radial enlargement of the 
dust aggregates after 400 and 3000 years. The x axis is 
the compact radii of the particles, the y axis is the ratio 
between the compact and the porous radius, this quantity 
is basically equal to \E'3. 



fluffy structures. However, when the distribution function 
relaxes in equilibrium, there are collisions between smaller 
and bigger aggregates as well which results in somewhat 
compacted aggregates. 



5.3.1. Results 

We performed a simulation with exactly the same initial 
conditions as in the last section but we included the poros- 
ity as an additional dust property in the model. The result 
can be seen in Figure [7] (the required CPU time here is 
also 1.5 hours). One can immediately see that including 
porosity increases the maximum particle mass by two or- 
ders of magnitude (five times larger particles in radius). 
This was al ready expected based on the work of Ormel 
et al. (|2007h . although due to rain out of bigger particles, 
they did not simulate particles bigger than 0.1 cm. 

We provide Figure [5] to give an impression how the 
porosity of the agglomerates change during the simulation. 
The x axis is the compact radius of the particles, the y axis 
is the ratio between the compact and the porous radii. 
This quantity is basically equal to \&3. Fractal growth is 
important for small particles creating fluffy agglomerates 
(until 10 -3 — 10~ 2 cm approximately), after this point 
the relative velocities become high enough so compact- 
ness becomes important. Before the particles reach a fully 
compacted stage they fragment, become monomers and a 
new cycle of growth starts. It is important to note that the 
porosity of the aggregates before the first fragmentation 
event is usually higher than the porosity values after equi- 
librium is reached. This can be seen in Figure [5] (grains 
after 400 years and 3000 years). The reason is that before 
the first fragmentation event, particles involved in colli- 
sions are typically equal sized so these particles produce 



5.3.2. Model comparison with Ormel et al. (2007) 



We compare our Mon te Carlo code with the one developed 
by Ormel et al. (|2007h . They use the Minimum Mass Solar 
Nebula disk model (MMSN) and somewhat different dust- 
parameters which we changed accordingly (distance from 
the central star = 1 AU, temperature = 280 K, density 
of the gas is 8.5 x 10~ 10 g/cm 3 , gas to dust ratio = 240, 
a = 10~ 4 ; monomer radius = 0.1 /im, monomer density 
= 3 g/cm 3 , surface energy density of the monomers = 25 
ergs/cm -2 ). 

They follow particle coagulation at one pressure scale 
height above the midplane of the disk. Because of this if 
the particles reach a critical stopping time {r ra i n = a/Q, 
where fi is the Kepler frequency), the particles rain out 
meaning that these particles leave the volume of the sim- 
ulation, the distribution function of the dust particles is 
collapsing as it can be se en in their figures (Figure 10 and 
11 in Ormel et al. (|2007l )V 

We do not include this effect in our model but we stop 
the simulation at the first rain out event and compare our 
distribution functions until this point. We use 10 4 repre- 
sentative particles (100 x 100) during the simulation. 

This can be seen at Figure El The reader is advised to 
exami ne thi s figure together with Figure 10. c. from Ormel 
et al. (|2007l ) because this is the figure we reproduced here. 
Furthermore we would like to point out that the scale of 
the y axis is different in the two figures. Our figure shows 
two orders of magnitude from the normalized distribution 



12 



Zsom & Dullemond: A representative particle approach 




Fig. 9. Distribution functions obtained by using Ormel et 
al. (|2007l ) input parameters. The continuous lines show 
the distribution functions at t = 10 years (thin line), 100 
years (thicker line), 1000 years (thickest line). The dotted 
line shows the distribution function at the time of the first 
compaction event (t = 1510 years), the dashed line shows 
the distribution function at the first rain out event (t = 
2900 years). 



functions whereas their figure covers more than 10 orders 
of magnitude from the real distribution function. 

Keeping these in mind, one can compare the results of 
the two Monte Carlo codes. 

The continuous lines at Figure [9] in this paper show 
the distribution functions at t = 10 years (thin line), 100 
years (thicker line), 1000 years (thickest line). The dotted 
line shows the distribution function at the time of the 
first compaction event (t = 1510 years), the dashed line 
shows the distribution function at the first rain out event 
(t = 2900 years). The same notation is used by Ormel et 
al. (|2007l ) at Figure 10. c. 

We compared the position of the peaks of the distribu- 
tion functions and the approximate shape of the curves. 
We can concl ude t hat our code reproduces the results of 
Ormel et al. (|2007h very well. 

The required CPU time to perform this simulation is 
only 10 minutes. One might ask why the CPU time is al- 
most ten times smaller now? Why do the previous simula- 
tions, which used the same number of representative par- 
ticles (10 4 ) and simulated approximately the same time 
interval (3000 years), take so long? The required CPU 
time does not scale linearly with the used number of par- 
ticles. It scales linearly with the number of collisions sim- 
ulated. The difference between this run and the previous 
two simulati ons is fragmentation. In the simulations of 
Ormel et al. (|2007l ) no fragmentation is happening because 
the growth timescales are longer. Using our initial param- 
eters, the first fragmentation event happens around 1000 
years, the number of small particles are never completely 
depleted after this time. As the small particles thereafter 



Fig. 10. The evolution of dust particles including the ef- 
fects of Brownian motion and turbulence, porosity and 
using two different monomer sizes (ai = 0.1/im and a-i = 
0.4/im). The particle distribution is saved after t — 3 x 10° 
years - dash-dot line, 3 x 10 1 years - dashed line, 3 x 10 2 
years - dotted line, and 3 x 10 3 years - continuous line. 



are always present, the number of collisions will be much 
higher than before. 

Also note that the porosities of these partic les would 
be smaller if the model of Ormel et al. (|2007| ) included 
fragmentation (for the reason see Sect. 15.3. Tj) . 



5.4. Monomer size distribution 

An interesting question which can easily be answered 
with our method is: How the mixture of different sized 
monomers change the maximum agglomerate size which 
can be reached? As we can see from Equations [14] and 
I15[ the rolling energy is lower for smaller monomers and 
of course the number of monomers in an agglomerate is 
much higher if the same agglomerate is built up of lighter 
monomers. This would mean higher fragmentation energy 
and one would expect that the particles would be harder 
to fragment resulting in bigger grains. 

We performed a simplified simulation to be able to 
answer this question. Only two different monomer sizes 
are considered here, a\ = O.l^tm and 02 = 0.4^m assuming 
that half of the mass (or representative particles) belongs 
to the small monomers, the other half belongs to the big 
monomers. 

One problem arises here with the rolling energy. The 
rolling energy changes with monomer size, and as our 
method cannot follow exactly the number of contacts in 
an aggregate and what kind of monomers are connected, 
we are forced to use an averaged rolling energy. One has to 
carefully consider what the average rolling energy should 
be. In our case, the big monomer is 64 times heavier 
than the small monomer. Let's assume that 50% of the 
mass of an aggregate is built up from small monomers; 
on the other hand, if we compare the number of differ- 
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ent monomers, the small monomers will be 64 times more 
numerous than the big ones. This means that the contri- 
bution of small monomers in the average rolling energy 
(Broil) should be higher. This can be achieved by using 
the following weighting: 

Eroll = r Erolh H r Eroll 1: (21) 

a\ + a 2 ai + a 2 

where E ro iii is the rolling energy between monomers with 
radius ai, E ro u 2 is the rolling energy between monomers 
with radius a 2 . 

As we can see in Figure 1 1 01 the maximum aggregate 
sizes reached are approximately an order of magnitude 
higher than on Figure [7] as it was predicted earlier in this 
Section. 

6. Conclusions and outlook 

We have shown that our representative particle method 
for aggregation of particles in astrophysical settings works 
well for standard kernels. It has the usual advantages of 
Monte Carlo methods that one can add particle properties 
easily and without loss of computational speed. Moreover, 
it naturally conserves the number of computational ele- 
ments, so there is no need to "add" or "remove" particles. 
Each representative particle represents a fixed portion of 
the total mass of solids. 

Our method may have various possible interesting ex- 
tensions and applications. Here we speculate on a few of 
these. For instance, the fact that each representative parti- 
cle corresponds to a fixed amount of solid mass makes the 
method ideal for implementation into spatially resolved 
models such as hydrodynamic simulations of planetary at- 
mospheres or protoplanetary disks. We can then follow 
the exact motion of each representative particle through 
the possibly turbulent environment, and thereby automat- 
ically treat the stochastic nature and deviation from a 
Boltzmann distribution of the motion of particles with 
stopping times of the same order as the turbulent eddy 
turn-over time. It is necessesary, however, to assure that 
a sufficiently large number of representative particles is 
present in each grid cell of the hydrodynamic simulation. 
For large scale hydrodynamic simulations this may lead 
to a very large computational demand for the coagulation 
computation, as well as for tracking the exact motion of 
these particles. If strong clumping of the particles hap- 
pens, however, much of the "action" anyway happens in 
these "clumps" , and it may then not be too critical that 
other grid cells are not sufficiently populated by represen- 
tative particles. This, however, is something that has to 
be experimented. 

Our representative particle method can in principle 
also be used to model the sublimation and condensation 
of dust grains. If a particle sublimates then the represen- 
tative particle becomes simply an atom or molecule of the 
vapor of this process. It will then follow the gas motion 
until the temperature becomes low enough that it can con- 
dense again. Other representative particles which are still 



in the solid phase may represent physical particles that 
can act as a condensation nucleus. Finally, in our method 
the properties of the particle can not only change due to 
collisions, but we can easily implement other environmen- 
tal factors in the alteration of particle properties. 

There are two main drawbacks of the method. First, it 
only works for large particle numbers, i.e. it cannot treat 
problems in which individual particles start dominating 
their immediate environment. Ormel's method and its ex- 
pected extension do not have this problem. Secondly, the 
method cannot be accelerated using implicit integration, 
while Brauer's method can. 

All in all we believe that this method may have in- 
teresting applications in the field of dust aggregation and 
droplet coagulation in protoplanetary disks and planetary 
atmospheres. 
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